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■ A Mathematics Notebook for Alpha Decay Using an Exponentially 
Diffuse Boundary Potential 



Eugene F. Chaffin and Daniel S. Banks 

Physics Department, Bob Jones University, 1700 Wade Hampton Blvd., 
Greenville, SC 29614 

Abstract Using the exponentially diffuse boundary potential of Green and Lee (1955), 
we develop a Mathematica notebook to treat alpha decay by using the complex energy 
approach of Preston (1947), as modified by Pierronne and Marquez (1978). Our 
approach allows the potential to change slowly from the well depth of the interior of the 
nucleus to the top of the Coulomb barrier, rather than the sudden step of the simple 
square well used by Pierronne and Marquez. Recent interest in possible variation in 
coupling constants such as the strong coupling a s motivates us to develop this 
algorithm which can allow numerical study of the variation of the decay constant of a 
nucleus such as U-238when the depth of the nuclear potential well changes. 

Introduction 

In recent years string theory the dependency of various constants on radii of 
compactified dimensions has made it appear that half-livesfor alpha decay may have 
been variable during the early history of the universe [1 -3,5,9-1 1 ]ln order to model this 
time dependence, we have used Mathematica to model the variation of the decay 
constant with change in depth of the potential well. More than 20 years ago, Pierronne 
and Marquez [6] treated the theory of alpha decay using a square well solution for the 
interior of the nucleus and coulombic solutions for the exterior, in a modification of the 
pioneer work by Preston [8] in which the complex nature of the alpha particle energy is 
utilized. The Pierronne and Marquez approach enables a finite well depth of around -30 
to -1 OOMeV to be used, whereas the earlier approach of Preston required an 
unrealistic assumption of a positive value for the well depth. The use of the square well 
leads to spherical Bessel functions for the interior solutions of the square well, which 
Pierronne and Marquez matched to repulsive coulomb solutions on the boundary of the 
potential well. 

In the 1950's Green and Lee [6] found solutions for a spherical well with an 
exponentially diffuse boundary potential V(r) = -V exp[(a-r)/a5],which were Bessel 
functions of nonintegral order for r>a, and the usual spherical Bessel functions for r<a, 
where a is the radius at which the exponential tail begins and 5 is a dimensionless 
parameter which characterizes the "shortness" of the tail. We have modified the 
approach of Pierronne and Marquez to use the Green and Lee solutions. Thus we 
match the logarithmic derivative of the spherical Bessel function solutions to the 
nonintegral order Bessel function solutions of Green and Lee at r = a, and then also 
match these solutions to the Coulombic wave functions at a larger radius r =b. 

The Pierronne and Marquez method also requires us to allow the one-bodyalpha 
particle energy to have a small complex part, thus modifying the bound states to allow 
tunneling through the barrier. Then we match the imaginary parts of the logarithmic 
derivative, which leads to the equations giving the decay constant. For the 
exponentially diffuse boundary wavefunctions of Green and Lee, the index v of the 
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Bessel function J v (kx) > becomes complex. In the limit where the imaginary part of the 
energy is small, a Taylor series expansion with v as the variable may be evaluated at 
the point where the imaginary part of v is zero, leading easily to the expression for the 
decay constant given in the Mathematica code below. 

Results 

The algorithm may be used to explore the variation of the decay constant as the well 
depth, nuclear radius, alpha particle energy at infinity, and other parameters vary. For 
example, Figure 1 below shows the negative logarithm of the decay constant plotted 
versus well depth, for a parent nucleus with Z=92, A= 238 (uranium-238).To produce 
this plot, the inner matching radius a was held constant, while the radius b at which the 
logarithmic derivatives of the Green and Lee solutions were matched to the coulombic 
wave functions was determined by Newton-Raphsoniteration, keeping the alpha 
particle energy fixed at 4.31 MeV. 



[Figure 1] 
Conclusions 

As has been clearly pointed out by Calmet and Fritsch [1 ,2], as the strong coupling constant as varies, several 
quantities may vary at once including the nucleon mass. Data such as the natural reactor at Oklo and the Sm-149 
cross section do not necessarily constrain these variations if more than one parameter varies at once. Due to the 
recent observations indicating the cosmological variation of the fine structure constant [10,1 1], this possibility has to 
be taken seriously. The algorithm given here may be used to explore the consequences of a variation of as on 
abundances of radioactive nuclides, such as those at Oklo. 
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■ The Mathematica Code 

(* 

An algorithm for finding the decay constant for variable 
well depths. The interior is modeled as a square well with 
an exponentially diffuse edge and a constant potential V-\ 
which matches the edge to the exterior Coulomb potential. 
Eugene F. Chaffin and Daniel S. Banks 
May 29, 2002 

The algorithm requires input giving the atomic number Z, 
the atomic mass number A,the kinetic energy of the emitted 
alpha particle E a and the depth V of the potential well. *) 
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(* The distance a is the point where the potential changes from the constant value 
for the interior of the square well to the exponentially diffuse potential Vi - 
VO Exp[(a-x)/(a <5) ] , in units of 10*10" 15 meters. 

The distance b is the point where the potential changes from the exponentially 
diffuse potential to the Coulomb potential, in units of 10*10~ 15 meters. *) 
a = 1; 
b = 1.045; 
VO = 50; 

V! = (8.98755* 10 9 * 2 * 90 * (1 . 6021 * 10" 19 ) 2 / (b*10~ 14 )) / (1 . 6021 * 10" 13 ) + 

V0 Exp[(a-b) / (a 5)]; 
6 = 0.035; 

M = 6.645 10 27 
Z = 90 
AA = 234 
E„ =4.31 

K ex = ( (2* M *E a *1.6021 10 13 ) 5 ) / (1.0546 10 ~ 34 ); 

Print ["K ex = ", K ex ] 
R = b*10" 14 

E = (1 . 05457 * 10" 34 ) 2 / (2 * 6 . 4 * 10" 27 * a 2 * 10" 28 ) / (1 . 6021 * 10" 13 ) 
/ VO 

60 = V it '• 

k = 2 * 5 * e ; 
RR = R / (10~ 14 ) ; 

v _ e ((l-RR/a)/(2.a)) . 

6.645 xlO" 27 
90 



234 
4.31 



K ex = 9.0836 xlO 14 
1.045 x 10" 14 
0.0542315 

This portion of the notebook provides an iterative process for finding the order 
n (non-integral)for the Bessel function 

which is the solution connecting the spherical Bessel function to the Coulomb 
function. The user provides an initial guess 
for n and calculates 

-2*5* (e ' / Tan [e' ] ) + 2 * 5 * e * BesselJ [n + 1, k]/BesselJ[n, k] . When 

the result is the same or 

close enough to the original value, the process is finished. 
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n = 1.60726 
e K = n / (2 * 6) 

e ' = -\] e 2 - e„ 2 + Vx 
k = 2 * S * e 

1.60726 
22.9609 
20.8182 
2.12548 

-2*6* (e ' / Tan[e ' ] ) + 2 * <5 * eo * BesselJ[n + 1, k] / BesselJ[n, k] 

1.60746 

77 = 0.6300*Z* ( (AA / (E a * (AA + 4) ) ) 5 ) ; 

GR[T_] := 

Exp [ - 2 * 77 * ((T* (1 - T) ) 5 + ArcSin[T° 5 ] - tt / 2) +0 .25*Log[T/ (1 - T) ] + 
(8*T 2 - 12 *T + 9) / (48 * Sqrt [T] * (1 - T) 1 5 ) / (2 *„) + 
(8*T - 3)/ (64*T*(1-T) 3 )/ ((2* r7 ) 2 )- 
(2048 *T 6 - 9216 *T 5 + 16128 *T 4 - 13440 *T 3 - 

12240 * T 2 + 7560 *T - 1890) / (92160 * T 1 5 * (1 - T) 45 ) / ( (2 *r]) 3 ) + 

3* (1024 *T 3 - 448*T 2 + 208 *T-39) / (8192* (T 2 ) * (1 - T) 6 ) / 

((2*,) 4 )- 

(- (262144 *T 10 - 1966080 *T 9 + 6389760 *T 8 - 11714560 * T 7 ) - 

(13178880 * T 6 - 9225216 * T 5 + 13520640 * T 4 ) - 

(-3588480 *T 3 + 2487240 * T 2 - 873180 * T + 130977) ) / (10321920 * (T 2 ' 5 ) * (1 - T) 7 5 ) / 

((2*n) 5 ) + ( (1105920 *T 5 - 55296 *T 4 + 314624 *T 3 - 159552 *T 2 ) ♦ 

(45576 *T - 5697) ) / (393216* (T 3 ) * (1 - T) 9 ) / ( (2 * r)) 6 ) ] ; 
GRD [T_] : = GR ' [ T ] ; 
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FR[T_] := 

.5 *Exp[ 2 * 77 * ( (T * (1 - T) ) 5 + ArcSin[T 5 ] - 7T / 2) + 
0.25*Log[T/ (1 - T) ] - 

((8*T 2 - 12 *T + 9) / (48*Sqrt[T] * (1-T) 1 ' 5 )) / 

(2*r7) + ( (8 * T - 3) / (64*T* (1-T) 3 )) / ( (2 * rj) 2 ) + 

( (2048 * T 6 - 9216 * T 5 + 16128 * T 4 - 13440 * T 3 - 12240 * T 2 + 

7560 *T- 1890) / (92160 *T 1 - 5 * (1-T) 45 )) / 

( (2 *r}) 3 ) + (3 * (1024 * T 3 - 448 * T 2 + 208 *T - 39) / 

(8192* (T 2 ) * (1 -T) 6 )) / ((2*?7) 4 ) + 

( (- (262144 *T 10 - 1966080 * T 9 + 6389760 * T 8 - 

11714560*T 7 ) - 

(13178880 *T 6 - 9225216 * T 5 + 13520640 * T 4 ) - 

(-3588480 *T 3 + 2487240 * T 2 - 873180 *T + 130977) ) / 

(10321920 * (T 2 - 5 ) * (1 - T) 7 - 5 ) ) / ( (2 *r]) 5 ) + 

( ( (1105920 * T 5 - 55296 * T 4 + 314624 * T 3 - 159552 * T 2 ) + 
(45576 *T - 5697) ) / 

(393216* (T 3 ) * (1-T) 9 )) / ((2*r7) 6 )]; 

FRD [T_] := FR ' [ T ] ; 

(* Having found n in the above procedure, 

we now start a loop for Newton-Raphson iteration to find the 
radius R by matching the logarithmic derivatives for the non- 
integral order Bessel function and the Coulomb wavefunction *) 
TES = 

i; 
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While TES > , 
L 1Q 20 

p e = K ex *R; 

RR = R / (10~ 14 ) ; 

v _ e ((l-RR/a)/(2»a)) . 

T = p e / (2* J?) ; 

(* When T is greater than 1, 

the Riccati method of expanding the Coulomb wavefunction is invalid 
(Froberg, 1955) , so an error message is printed in that case *) 
If[T>l, Print ["T = ", T] ] ; 

FOP = FRD [T] / (2 * 77) ; 
GOP = GRD [T] / (2 * 77) ; 

(* Find the Wronskian as a check on the accuracy *) 

W = GR[T] *FOP - FR[T] *GOP; 

If [Abs[l - W] > 10" 4 , Print ["W = " , W] ] ; 

GOPP = - GR [ T ] * (1 - rj *2 /p e ) ; 

FN = K ex *GOP / GR [T] + 

1 / (2 * 5 *a* 1(T 14 ) * (n - k* v*BesselJ[n + 1, k* v] / BesselJ[n, k* v] ) ; 

FNP = (K ex * GOPP / GR [T] - K ex * GOP * GOP / (GR [T] 2 ) ) * K ex + 

(1 / (2 *5*a*10" 14 )) * 
( (k 2 * v 2 / (2 * 6 *a* 10~ 14 *BesselJ[n, k* v] ) ) * ( (2 / (k* v) ) *BesselJ[n + 1, k * v] + 
(BesselJ[n + 1, k * v] ) 2 / BesselJ[n, k * v] - BesselJfn + 2, k * v] ) ) ; 

RP = R - FN /FNP; 
TES = Abs [RP - R] ; 

R = RP;] 

(* end of Newton-Raphson iteration loop *) 
Print ["R = " , R] 

R= 1.04499 x 10" 14 

B = (K ex *RP* ( GR [ T ] * GOPP -GOP 2 ) + GR [ T ] *GOP) / (GR [T] 2 ) ; 
CC = B* ((2*M/ (E a * 1.6021 10" 13 )) ') / (1.0546 10" 34 ) / A; 

(1-R/(,.10-")) 

x = k * e 2.6 ■ 



v = n; 
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Q = e 



(l-R/ (a*10~ 14 ) ) /2*S 



1/X + 



— *BesselJ[v, x] - BesselJ[v + l, x] 



(BesselJ[v, x] ) ' 
2" v *PolyGamma[v + 1] * x v 2" v " 2 * PolyGamma [v + 2] * x v+2 



x*BesselJ[v, x] 









Gamma [v + 1] 




Gamma [ v + 2 ] 


2 






* PolyGamma [v + 3] * 


x v+4 


2 _v_6 * PolyGamma [v + 4] * x v+6 
+ 








Gamma [ v + 3 ] 




Gamma [v + 4] 


2 




g 


* PolyGamma [v + 5] * 


v+fl 
X V + S 


2 -v-io # PolyGamma[v + 6] * x v+1 ° 
+ 








Gamma [ v + 5 ] 




Gamma [v + 6] 


2 


-V- 


12 


* PolyGamma [v + 7] * 


x v + 12 


2" v - 14 * PolyGamma [v + 8] * x v+14 
- + 








Gamma [ v + 7 ] 




Gamma [v + 8] 


2 


-V- 


16 


* PolyGamma [ v + 9 ] * 


x v + 16 


2" v - 18 *PolyGamma[v + 10] * x v+18 
- + 








Gamma [ v + 9 ] 




Gamma [v + 10] 


2 


-V- 


20 


*PolyGamma[v + 11] 


* X v+2 ° 


2 -v-22 * p iy Gamma [ V + 12] * x v+22 








Gamma[v + 11] 




Gamma [v + 12] 


2 


-V- 


24 


* PolyGamma [v + 13] 


* X v+24 


2" v - 26 * PolyGamma [v + 14] * x v+26 



Gamma [v + 13] 
2" v " 28 * PolyGamma [v + 15] * x v+2S )' 
Gamma [v + 15] 

(* Calculate the decay constant *) 



Gamma [v + 14] 



X = 



CC + 



k*<5*M*a*10 



*Q 



* (1.0546 10" 34 ) 



(1.0546 10" 34 ) 
( GR [ T ] * GR [ T ] ) ; 

Print [ "A = " , A] 

A= 1.6956x10"" 
Vi 

38.6247 

(8.98755* 10 9 *2 * 90 * (1 . 6021 * 10" 19 ) 2 / (b* 10 -14 )) / (1 . 6021 * 10" 13 ) 

24.802 

(* Print the radius resulting from the matching *) 
Print ["R = " , R] 

R= 1.04499 x 10" 14 
Q 

1.12317 

cc 



-2.78579 x 1(T 
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k 

2.12548 
M 

6.645 x 10" 27 
6 

0.035 

a 

1 

6 ' = -\j 6 2 - e„ 2 + Vi 

20.8182 

h[p_] := ( (p) 05 ) BesselJ[0.5, e' *p]; 
h[l] 

0.161212 

BesselJ[n, k] 

0.482814 

h[l] /BesselJ[n, k] 

0.333902 

| (((p) 5 ) BesselJ[0.5, e' *p]) 2 (2p 
Jo 

0.0155523 

f 10 2 
N[ (0.334227) 2 (BesselJ[n, k*E (1 " p)/(2 * 6) ] ) dp] 

0.000842111 

1 /VO. 015553 + .000842111 

7.80985 

7.80985* (0.334227) 

2.61026 
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GR[T_] := (8.51346* 10" 12 ) (2/ (a * 10~ 14 ) ) °' 5 * ((2*r)*T) /K ex ) * 

Exp [ - 2 * r) * ((T* (1 - T) ) 5 + ArcSin[T 5 ] - tt/2) +o .25 *Log[T/ (1 - T) ] + 
(8*T 2 - 12 *T + 9) / (48*Sqrt [T] * (1 - T) 1 5 ) / (2* n ) + 
(8*T - 3)/ (64*T*(1-T) 3 )/ ( (2 *r?) 2 ) - 
(2048 *T 6 - 9216 *T S + 16128 *T 4 - 13440 *T 3 - 

12240 * T 2 + 7560 *T - 1890) / (92160 * T 1 5 * (1 - T) 4 5 ) / ( (2 * rj) 3 ) + 

3* (1024 *T 3 - 448* T 2 + 208 *T-39) / (8192* (T 2 ) * (1 - T) 6 ) / 
((2*,) 4 )- 

(- (262144 *T 10 - 1966080 *T 9 + 6389760 *T 8 - 11714560 *T 7 ) - 

(13178880 * T 6 - 9225216 * T 5 + 13520640 * T 4 ) - 

(-3588480 *T 3 + 2487240 * T 2 - 873180 * T + 130977) ) / (10321920* (T 2 S ) * (1 - T) 7 5 ) / 

( (2*77) 5 ) + ( (1105920 *T 5 -55296*T 4 +314624*T 3 -159552*T 2 ) + 

(45576 *T - 5697) ) / (393216* (T 3 ) * (1 - T) 9 ) / ( (2 * rj) 6 ) ] 
TR = R*K ex / (2*!7) 

0.175257 

GR [ TR] 

0.603691 

gt[p_] := 2.61006*BesselJ[n, k * E (1 " p) ' (2 * i) ] / p; 
pR = R/ 10" 14 

General:: spell 1 : Possible spelling error: new symbol name "pR" is similar to existing symbol "p". 

1.04499 

gt[pR] 

0.603773 

0.603888/ (7.09333* 10 10 ) 

8.51346 xlO" 12 



AlphaDiffuse. nb 



11 



g[p_] := Which[p < 1, (7.80923) * ((p) 05 ) BesselJ[0. 5, e' *p], 1 <, p < (R/10 -14 ), 

2.61006*BesselJ[n, k * E (1 " p) ' i2 * s) ] I p, (R / (a*10~ 14 )) < p, GR[K ex p*a*10 -14 / (2*77)]]; 
fV[p_] :=Which[p<l, -VO, 

l<p<b, V! - VO Exp [ (1 - p) / ( <5) ] , p > b, 

(8.98755* 10 9 *2 * 90 * (1 . 6021 * 10" 19 ) 2 / (p * 10" 14 ) ) / (1.6021 * 10" 13 ) ] ; 

Plot[{fV[p], 10*g[p]}, {p, 0.01, 2}, PlotRange -> All, PlotLabel -> "\ ! \ (\* 
StyleBox[\"Wavefunction\", \nFontFamily->\ "Arial\ " , \nFontSize->24] \) \ !\ (\* 
StyleBox [ \ " x\ " , \nFontFamily- >\ "Arial\ " , \nFontSize- >2 4 ] \) \ ! \ ( \* 
Sty leBox [ \ " 1 \ " , \nFontFamily- > \ " Arial\ " , \nFont Si ze- >24]\)\!\(\* 
StyleBox [ \ " \ " , \nFontFamily- >\ " Arial\ " , \nFontSize- >2 4 ] \) \ ! \ ( \ * 
StyleBox [V'andV , \nFontFamily->\ "Arial\ " , \nFontSize->24] \) \ ! \ (\* 
StyleBox [ \ " \ " , \nFontFamily- >\ " Arial\ " , \nFontSize- >2 4 ] \) \ ! \ ( \ * 
StyleBox [\"Potential\" , \nFontFamily->\"Arial\" , \nFontSize->24] \) \ ! \ (\* 
StyleBox [ \ " \ " , \nFontFamily- >\ " Arial\ " , \nFontSize- >2 4 ] \) \ ! \ ( \ * 
StyleBox [\"versus\" , \nFontFamily->\"Arial\" , \nFontSize->24] \) \ ! \ (\* 
StyleBox [ \ " \ " , \nFontFamily- >\ " Arial\ " , \nFontSize- >2 4 ] \) \ ! \ ( \ * 
StyleBox [\ "Distance\ ", \nFontFamily->\ "Arial\ " , \nFontSize->24] \) ", 

AxesLabel -> {p , psi[V]}, PlotStyle ^ {RGBColor [1, 0, 0], RGBColorfO, 1, 0]}, 

TextStyle -» {FontSize -» 24} ] ; 
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